*
*
* GARCH MC for delta-method (all 4 experiments)
* 5/20/24
*	
* -----------------------------------------------------------------------
cd "~/Dropbox/Benton-Jordan-Philips/final-JOP-replication"

* run garch sim program:
do "scripts-programs/synthgarch.ado"

* -----------------------------------------------------------------------------
* -----------------------------------------------------------------------------
* 					---------- SCENARIO 1 ----------------------
* -----------------------------------------------------------------------------
/*	ARCH = 0.1
	GARCH = 0.5
	T = 250
	beta_x = .5
	beta_z = -1
	constant = 0.5
	phi = 0.25
	x is continuous, z is dichotomous, epsilon is normal
	burnin = 250
	sims = 100
	scenname = scenario1
*/

loc dgp_sims = 500 // how many draws from the DGP are we taking?
loc rep_sims = 500 // how many reps OF a dgp_sim do we want to take?

* -----------------------------------------------------------------------
* -->	STEP 1: GENERATE SERIES. 
* this is done for all steps below
synthgarch, time(250) garchparam(0.5) ///
 archparam(0.1) xparam(.5) zparam(-1) consparam(0.5) ///
 phiparam(0.25) burnin(250) seedz(2000221) sims(`dgp_sims') scenname(scenario1)
* -----------------------------------------------------------------------

* * -----------------------------------------------------------------------
* -->	STEP 2: DELTA METHOD
* We're producing 1 prediction per sim. 
* estimate GARCH model on dataset m. 
* Use delta method to create sigma2 prediction for m. Burn in to make stable prediction
* record prediction, SE, UL and LL
mat storage = J(`dgp_sims',4,.) // cols are:
* 1. delta method prediction
* 2. ll prediction
* 3. ul prediction
* 4. sd prediction

use "dgp_scenario1.dta", clear
set seed 1245039932
tsset t
set obs 2000 // ensure we can store up to M results

* now estimate M GARCH models, saving coefficients from HET eq:
forv i = 1/`dgp_sims' { 
	* run garch
	cap arch y_`i' l.y_`i' x_`i' z_`i', arch(1) garch(1) het(x_`i' z_`i')
	
	if "`e(converged)'" == "1" {
		* init guess, then burn-in to create stable prediction:
		cap noisily nlcom exp(_b[HET:_cons] + _b[HET:x_`i']*1 + _b[HET:z_`i']*1)
		if _rc == 0 { // only continue if nlcom success
			forv z = 1/100 {
				mat delta_b = r(b)
				local delta_b = delta_b[1,1]
				cap nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
					_b[HET:x_`i']*1 + _b[HET:z_`i']*1) + _b[ARCH:L.arch]*0
				* last loop, save values (assuming nlcom success)
				if _rc == 0 & `z' == 100 { 
					* now that we're at a stable expected sigma2, obtain UL,
					* LL and the prediction
					mat delta_fin = r(table)
					mat storage[`i',1] = delta_fin[1,1] // expected value
					mat storage[`i',2] = delta_fin[2,1] // SD
					mat storage[`i',3] = delta_fin[5,1] // LL
					mat storage[`i',4] = delta_fin[6,1] // UL
				}
			}
			
			
		}	
	}
	else {
		* do nothing. no convergence
	}
}
mat list storage
clear
svmat storage
* 1. estimated prediction
rename storage1 pred_DELTA
* 2. estimated SD
rename storage2 pred_sd_DELTA
* 3. estimated LL
rename storage3 pred_ll_DELTA
* 4. estimated UL
rename storage4 pred_ul_DELTA
gen sims = _n
compress
save "data/dgp_scenario1_DELTAPRED.dta", replace
* -----------------------------------------------------------------------


* -----------------------------------------------------------------------------
* -----------------------------------------------------------------------------
* ---------- SCENARIO 2 (T=1000, normal dist resids) ----------------------
* -----------------------------------------------------------------------------
/*	ARCH = 0.1
	GARCH = 0.5
	T = 1000
	beta_x = .5
	beta_z = -1
	constant = 0.5
	phi = 0.25
	x is continuous, z is dichotomous, epsilon is normal
	burnin = 250
	sims = 100
	scenname = scenario2
*/

loc dgp_sims = 500 // how many draws from the DGP are we taking?
loc rep_sims = 500 // how many reps OF a dgp_sim do we want to take?

* -----------------------------------------------------------------------
* -->	STEP 1: GENERATE SERIES. 
* this is done for all steps below
synthgarch, time(1000) garchparam(0.5) ///
 archparam(0.1) xparam(.5) zparam(-1) consparam(0.5) ///
 phiparam(0.25) burnin(250) seedz(43055029) sims(`dgp_sims') ///
 scenname(scenario2)
* -----------------------------------------------------------------------

* * -----------------------------------------------------------------------
* -->	STEP 2: DELTA METHOD
* We're producing 1 prediction per sim. 
* estimate GARCH model on dataset m. 
* Use delta method to create sigma2 prediction for m. Burn in to make stable prediction
* record prediction, SE, UL and LL
mat storage = J(`dgp_sims',4,.) // cols are:
* 1. delta method prediction
* 2. ll prediction
* 3. ul prediction
* 4. sd prediction

use "dgp_scenario2.dta", clear
set seed 409090
tsset t
set obs 2000 // ensure we can store up to M results

* now estimate M GARCH models, saving coefficients from HET eq:
forv i = 1/`dgp_sims' { 
	* run garch
	cap arch y_`i' l.y_`i' x_`i' z_`i', arch(1) garch(1) het(x_`i' z_`i')
	
	if "`e(converged)'" == "1" {
		* init guess, then burn-in to create stable prediction:
		cap noisily nlcom exp(_b[HET:_cons] + _b[HET:x_`i']*1 + _b[HET:z_`i']*1)
		if _rc == 0 { // only continue if nlcom success
			forv z = 1/100 {
				mat delta_b = r(b)
				local delta_b = delta_b[1,1]
				cap nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
					_b[HET:x_`i']*1 + _b[HET:z_`i']*1) + _b[ARCH:L.arch]*0
				* last loop, save values (assuming nlcom success)
				if _rc == 0 & `z' == 100 { 
					* now that we're at a stable expected sigma2, obtain UL,
					* LL and the prediction
					mat delta_fin = r(table)
					mat storage[`i',1] = delta_fin[1,1] // expected value
					mat storage[`i',2] = delta_fin[2,1] // SD
					mat storage[`i',3] = delta_fin[5,1] // LL
					mat storage[`i',4] = delta_fin[6,1] // UL
				}
			}
			
			
		}	
	}
	else {
		* do nothing. no convergence
	}
}
mat list storage
clear
svmat storage
* 1. estimated prediction
rename storage1 pred_DELTA
* 2. estimated SD
rename storage2 pred_sd_DELTA
* 3. estimated LL
rename storage3 pred_ll_DELTA
* 4. estimated UL
rename storage4 pred_ul_DELTA
gen sims = _n
compress
save "data/dgp_scenario2_DELTAPRED.dta", replace
* -----------------------------------------------------------------------






* -----------------------------------------------------------------------------
* -----------------------------------------------------------------------------
* 					---------- SCENARIO 1-t ----------------------
* -----------------------------------------------------------------------------
/*	ARCH = 0.1
	GARCH = 0.5
	T = 250
	beta_x = .5
	beta_z = -1
	constant = 0.5
	phi = 0.25
	x is continuous, z is dichotomous, epsilon is normal
	burnin = 250
	sims = 100
	scenname = scenario1t
*/

loc dgp_sims = 500 // how many draws from the DGP are we taking?
loc rep_sims = 500 // how many reps OF a dgp_sim do we want to take?

* -----------------------------------------------------------------------
* -->	STEP 1: GENERATE SERIES. 
* this is done for all steps below
synthgarch, time(250) garchparam(0.5) ///
 archparam(0.1) xparam(.5) zparam(-1) consparam(0.5) ///
 phiparam(0.25) burnin(250) seedz(430029) sims(`dgp_sims') ///
 scenname(scenario1t) tdist
* -----------------------------------------------------------------------

* * -----------------------------------------------------------------------
* -->	STEP 2: DELTA METHOD
* We're producing 1 prediction per sim. 
* estimate GARCH model on dataset m. 
* Use delta method to create sigma2 prediction for m. Burn in to make stable prediction
* record prediction, SE, UL and LL
mat storage = J(`dgp_sims',4,.) // cols are:
* 1. delta method prediction
* 2. ll prediction
* 3. ul prediction
* 4. sd prediction

use "dgp_scenario1t.dta", clear
set seed 4898999
tsset t
set obs 2000 // ensure we can store up to M results

* now estimate M GARCH models, saving coefficients from HET eq:
forv i = 1/`dgp_sims' { 
	* run garch
	cap arch y_`i' l.y_`i' x_`i' z_`i', arch(1) garch(1) het(x_`i' z_`i')
	
	if "`e(converged)'" == "1" {
		* init guess, then burn-in to create stable prediction:
		cap noisily nlcom exp(_b[HET:_cons] + _b[HET:x_`i']*1 + _b[HET:z_`i']*1)
		if _rc == 0 { // only continue if nlcom success
			forv z = 1/100 {
				mat delta_b = r(b)
				local delta_b = delta_b[1,1]
				cap nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
					_b[HET:x_`i']*1 + _b[HET:z_`i']*1) + _b[ARCH:L.arch]*0
				* last loop, save values (assuming nlcom success)
				if _rc == 0 & `z' == 100 { 
					* now that we're at a stable expected sigma2, obtain UL,
					* LL and the prediction
					mat delta_fin = r(table)
					mat storage[`i',1] = delta_fin[1,1] // expected value
					mat storage[`i',2] = delta_fin[2,1] // SD
					mat storage[`i',3] = delta_fin[5,1] // LL
					mat storage[`i',4] = delta_fin[6,1] // UL
				}
			}
			
			
		}	
	}
	else {
		* do nothing. no convergence
	}
}
mat list storage
clear
svmat storage
* 1. estimated prediction
rename storage1 pred_DELTA
* 2. estimated SD
rename storage2 pred_sd_DELTA
* 3. estimated LL
rename storage3 pred_ll_DELTA
* 4. estimated UL
rename storage4 pred_ul_DELTA
gen sims = _n
compress
save "data/dgp_scenario1t_DELTAPRED.dta", replace
* -----------------------------------------------------------------------




* -----------------------------------------------------------------------------
* -----------------------------------------------------------------------------
* 					---------- SCENARIO 2 (T=1000, t dist resids) ----------------------
* -----------------------------------------------------------------------------
/*	ARCH = 0.1
	GARCH = 0.5
	T = 1000
	beta_x = .5
	beta_z = -1
	constant = 0.5
	phi = 0.25
	x is continuous, z is dichotomous, epsilon is normal
	burnin = 250
	sims = 100
	scenname = scenario2t
*/

loc dgp_sims = 500 // how many draws from the DGP are we taking?
loc rep_sims = 500 // how many reps OF a dgp_sim do we want to take?

* -----------------------------------------------------------------------
* -->	STEP 1: GENERATE SERIES. 
* this is done for all steps below
synthgarch, time(1000) garchparam(0.5) ///
 archparam(0.1) xparam(.5) zparam(-1) consparam(0.5) ///
 phiparam(0.25) burnin(250) seedz(43055029) sims(`dgp_sims') ///
 scenname(scenario2t) tdist
* -----------------------------------------------------------------------

* * -----------------------------------------------------------------------
* -->	STEP 2: DELTA METHOD
* We're producing 1 prediction per sim. 
* estimate GARCH model on dataset m. 
* Use delta method to create sigma2 prediction for m. Burn in to make stable prediction
* record prediction, SE, UL and LL
mat storage = J(`dgp_sims',4,.) // cols are:
* 1. delta method prediction
* 2. ll prediction
* 3. ul prediction
* 4. sd prediction

use "dgp_scenario2t.dta", clear
set seed 4353454
tsset t
set obs 2000 // ensure we can store up to M results

* now estimate M GARCH models, saving coefficients from HET eq:
forv i = 1/`dgp_sims' { 
	* run garch
	cap arch y_`i' l.y_`i' x_`i' z_`i', arch(1) garch(1) het(x_`i' z_`i')
	
	if "`e(converged)'" == "1" {
		* init guess, then burn-in to create stable prediction:
		cap noisily nlcom exp(_b[HET:_cons] + _b[HET:x_`i']*1 + _b[HET:z_`i']*1)
		if _rc == 0 { // only continue if nlcom success
			forv z = 1/100 {
				mat delta_b = r(b)
				local delta_b = delta_b[1,1]
				cap nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
					_b[HET:x_`i']*1 + _b[HET:z_`i']*1) + _b[ARCH:L.arch]*0
				* last loop, save values (assuming nlcom success)
				if _rc == 0 & `z' == 100 { 
					* now that we're at a stable expected sigma2, obtain UL,
					* LL and the prediction
					mat delta_fin = r(table)
					mat storage[`i',1] = delta_fin[1,1] // expected value
					mat storage[`i',2] = delta_fin[2,1] // SD
					mat storage[`i',3] = delta_fin[5,1] // LL
					mat storage[`i',4] = delta_fin[6,1] // UL
				}
			}
			
			
		}	
	}
	else {
		* do nothing. no convergence
	}
}
mat list storage
clear
svmat storage
* 1. estimated prediction
rename storage1 pred_DELTA
* 2. estimated SD
rename storage2 pred_sd_DELTA
* 3. estimated LL
rename storage3 pred_ll_DELTA
* 4. estimated UL
rename storage4 pred_ul_DELTA
gen sims = _n
compress
save "data/dgp_scenario2t_DELTAPRED.dta", replace
* -----------------------------------------------------------------------
